function [A, b] = process_boundary_2d(A, b, mesh, pde)
%% Dirichlet
for k = 1:size(pde.bdary,2)
    if 0 == pde.bdary(1,k)
        bdn = boundarynodes(mesh, k);
        A(bdn,:) = 0;
        A(sub2ind(size(A), bdn, bdn)) = 1;
        b(bdn) = pde.u1(mesh.node(bdn,1),mesh.node(bdn,2));
        bdn = bdn + mesh.Nb;
        A(bdn,:) = 0;
        A(sub2ind(size(A), bdn, bdn)) = 1;
        b(bdn) = pde.u2(mesh.node(bdn-mesh.Nb,1),mesh.node(bdn-mesh.Nb,2));
    end
end
end

%% boundarynodes
function bdn = boundarynodes(mesh, index)
switch index
    case 1
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(1:2*mesh.Nx+1);
        else
            bdn = mesh.bdn(1:mesh.Nx+1);
        end
    case 2
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(2*mesh.Nx+1:2*mesh.Nx+2*mesh.Ny+1);
        else
            bdn = mesh.bdn(mesh.Nx+1:mesh.Nx+1+mesh.Ny);
        end
    case 3
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(2*mesh.Nx+2*mesh.Ny+1:4*mesh.Nx+2*mesh.Ny+1);
        else
            bdn = mesh.bdn(mesh.Nx+1+mesh.Ny:2*mesh.Nx+1+mesh.Ny);
        end
    case 4
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(4*mesh.Nx+2*mesh.Ny+1:end);
        else
            bdn = mesh.bdn(2*mesh.Nx+1+mesh.Ny:end);
        end
        bdn = [bdn; 1];
    otherwise
        error("Invalid edge number");
end
end